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Nomenclature  (cont.) 


_  K  .  t  i « 

-u'v'  Reynolds  stress  in  the  boundary  layer 

P 


x  -  arc  length  distance  along  body  =  — 

L 

* 

y  -  distance  normal  to  body  surface  =  ^ 

L 

*  * 

£  -  turbulent  eddy  viscosity  =  -  ■■ 

U 

<()  -  body  surface  angle,  see  Figure  2 

* 

^  -  Stokes  stream  function  =  — £-~~r — j- 

p  L  U 


ip  -  Stokes  stream  function  at  the  boundary-layer  edge  = 
o 


*  *2  * 
p  L  U 


T  -  shear  stress  = 


*  *2 

P  Uoo 


All  of  the  above  definitions  are  nondimensional  and  based  on  the  following 
dimensional  reference  quantities. 


-  reference  length 

-  dynamic  viscosity 

-  density 


-  shear  stress 


free  stream  velocity 
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INTRODUCTION 


1.1  Origin  and  Importance  of  the  Investigation 

In  dealing  with  the  flow  field  around  axisymmetric  bodies,  difficulties 
arise  in  the  vicinity  of  the  tail.  Here  the  thickness  of  the  boundary  layer 
produces  a  displacement  effect,  altering  the  inviscid  pressure  distribution. 
The  resultant  curvature  of  the  outer  streamlines  causes  the  development  of 
a  significant  normal  pressure  gradient  across  the  boundary  layer,  thus 
negating  the  thin  layer  approximation  currently  used  to  calculate  boundary 
layer  flow.  This  is  the  so-called  strong-interaction  region,  named  for 
the  strongly  coupled  turbulent  and  inviscid  flow  that  exists  at  the  aft 
end  of  the  body.  The  accuracy  of  the  prediction  of  the  flow  field  in 
this  area  is  important  due  to  its  major  contribution  to  the  total  flow  drag. 

Present  treatment  of  the  strong-interaction  phenomenon  involves  the 
simultaneous  solution  of  the  equations  governing  the  boundary  layer  and 
inviscid  flow  fields  using  a  displacement  body  technique.  The  procedure 
is  iterative.  The  viscous  flow  is  calculated  using  a  pressure  distribution 
determined  from  inviscid  flow  theory.  The  displacement  thickness  is 
calculated  and  added  to  the  original  body.  This  displacement  body  is  then 
used  to  calculate  a  new  inviscid  pressure  distribution  and  the  cycle  is 
repeated  until  convergence.  An  example  of  this  method  is  presented  by 
Hoffman  [1]*. 


1.2  Statement  of  the  Problem 


In  the  correct  application  of  the  displacement  body  method  the  normal 
pressure  gradient  is  accounted  for  by  mapping  the  pressure  distribution  of 
the  displacement  body  back  to  the  original  body  while  conserving  the  arc 
length.  Whether  or  not  this  properly  models  the  strong-interaction  flow 
has  become  a  point  of  contention  [2],  By  directly  including  the  normal 
pressure  gradient  the  accuracy  of  the  flow  field  prediction  should  be 
improved . 

The  mechanics  of  this  problem  requires  the  coupling  of  the  normal 
momentum  equation  to  the  streamwise  momentum  equation  and  the  continuity 
equation.  All  simplifications  with  regard  to  the  pressure  terms  in 
these  equations  are  dropped.  As  in  the  displacement  body  technique,  the 
simultaneous  solution  of  viscous  and  inviscid  equations  involves  an 
iterative  procedure,  but  there  are  two  major  differences.  First,  the 
potential  flow  solution  is  patched  to  the  edge  of  the  boundary  layer 
rather  than  the  displacement  body.  Second,  the  pressure  variation  through 
the  boundary  layer  is  calculated  at  each  station  in  the  flow  using  the 
edge  pressure  determined  from  the  potential  flow. 


Numbers  in  brackets  designate  References  at  the  end  of  report. 
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The  final  results  will  be  the  prediction  of  the  pressure  distributions 
on  the  body,  at  the  edge  of  the  boundary  layer  and  wake,  and  along  the  wake 
centerline.  Once  these  values  are  known,  conventional  methods  can  be 
employed  to  calculate  the  important  flow  parameters,  including  the  drag 
coefficient  of  the  body. 


1.3  General  Scope  of  the  Investigation 

The  subject  of  the  investigation  is  the  F-57  Body  examined  by  Patel 
and  Lee  [3],  shown  in  Figure  1.  This  body  was  chosen  for  several  reasons,  most 
important  of  which  is  its  cusped  tail.  The  cusp  allows  the  boundary-layer 
calculations  to  continue  smoothly  into  the  wake  without  the  discontinuity 
that  would  accompany  a  finite-angle  tail.  Also,  the  large  amount  of 
experimental  data  already  available  on  the  F-57  Body  makes  it  an  excellent 
subject  for  the  initial  testing  of  a  numerical  procedure.  After  obtaining 
a  converged  solution  for  the  simplified  case,  the  method  could  be  extended 
to  include  finite-angle  tails.  Due  to  insufficient  time,  that  problem  was 
not  addressed  in  this  report. 


THEORETICAL  DEVELOPMENT 


2.1  Boundary  Layer  Theory 


The  boundary  layer  equations  for  a  body  of  revolution,  considering 
only  incompressible,  steady  flow,  are  [4]: 


Continuity: 


X-Momentum : 


Y-Momentum: 


h (ru)  +  h (rv)  =  0  » 


3u  3u  3p  .  1  3 

U  -  +  v  ■  -  =  —  4-  —  - 

3x  3y  3x  r  3y 


u 


3v  3\r  _  _  3£ 
3x  V  3y  3y 


> 


■ 

■  ' 

r 

1  3y  t — i 

— —  ~~  ~  U  V 

Re  3y 

• 

(1) 

(2) 

(3) 


where  all  variables  are  dimensionless.  The  coordinate  system  employed  is 
shown  in  Figure  2  as  it  applies  to  the  aft  portion  of  the  body.  From 
Figure  2,  the  relationship  between  r  and  y  is 


r(x,y)  =  rQ(x)  +  y  cos  <fi  (4) 

Equations  (1)  -  (3)  include  the  assumption  that  all  the  stress  terms  in 
the  original  Navier-Stokes  equations,  except  the  aXy  gradient,  are 
negligibly  small.  They  do  not,  however,  allow  for  the  pressure  gradient 
normal  to  the  surface  to  go  to  zero.  The  boundary  conditions  are: 


u(x,0)  =  0 
v(x,0)  ■  0 


lim  u(x,y)  -  u  (x) 
e 

y-H» 


(5) 

(6) 
(7) 
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The  Reynolds  stress  is  related  to  the  velocity  field  by  assuming  the 
relationship 

— i — T  e  3u  /Q. 

”u  v  =  Re  3y  (8) 

where  e  is  the  eddy  viscosity.  Using  Equation  (7),  the  x-momentum  equation 
becomes 

3u  ,  3u  3p  ,  1  3  [/,  ,  .  3u  1 

“  a?  +  v  a?’ '  + TSa?  lr(l  +  E)  sfj  '  (9) 

At  this  point,  the  equations  must  be  put  into  a  form  compatible  with  the 
existing  ARL  boundary  layer  program  in  scaled  physical  variables,  which 
employs  Keller’s  Box  Method.  The  Reynolds  number  is  eliminated  from  the 
streamwise  momentum  equation  by  introducing  the  following  scaled  variables: 


•  -X- 


Equations  (1),  (2),  and  (3)  respectively 


h  +  h  (tv)  ■  0 


3u  .  3u  3p  .13 

u  —  +  v  —  =  -r*-  + - 

3x  —  3x  r  — 

3y  3y 


r (1  +  e) 


3v  .  3v  D 

u  —  +  v  —  =  -  Re 

3x  3y 


i£ 

_  9 


In  order  to  use  the  box  method,  the  Stokes  stream  function  is  introduced 
to  manipulate  the  equations  into  the  necessary  first-order  form.  The 
scaled  Stokes  stream  function  is  defined  as 


rv  =  - 

-  ’  rv  3x 

3y 
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Equation  (12)  is  Identically  satisfied  by  the  stream  function.  Grouping 
terms  and  defining  the  laminar  shear  stress  as 


T 


3u 

3y 


and 


b  =  (1  +  e) 


(16) 


(17) 


the  streamwise  momentum  equation.  Equation  (13),  can  be  written  in  first- 
order  form. 


f  h  (u2  +  2p)  ■  H T  ■  7= (rbT)  = 

3y 


(18) 


The  first-order  set  of  equations  to  be  solved  is: 


3  /  .  x  r  3  ,  2  ,  „  . 

—  (rbr)  -  Y  (u  +  2p)  -  ^  t 


3y 

3y 

jht 

ay 


3x 


=  T 


ru 


(19) 


(20) 


(21) 


,-2, 


1  3 (v  )  ,  D  3p  3v 

2  37  Re  sy  ~  "  u  3x 

The  boundary  conditions  for  this  set  of  equations  are: 


(22) 


u  =  v  =  ^  =  0 

at 

y  = 

0 

(23) 

U  “►  u 

e 

at 

7  -*• 

do 

(24) 

T  *  0 

at 

7  - 

0  (only  in  wake) 

(25) 

mmm 


-10- 


24  July  1980 
MSD:cag 


Equations  (19)  -  (21)  form  the  conventional,  parabolic  set  of  boundary-layer 
equations.  The  addition  of  Equation  (22),  however,  couples  this  set  to  the 
elliptic,  potential  flow  equations.  Fortunately,  this  coupling  is  weak, 
thus  allowing  the  normal  momentum  equation  to  lag  one  cycle  in  an  iterative 
calculation  scheme. 


2.2  Eddy  Viscosity  Model 

The  algebraic  eddy  viscosity  model  used  in  the  boundary-layer  calcula¬ 
tions  is  a  modification  of  the  Baldwin-Lomax  turbulence  model  [6],  [7]. 

It  consists  of  an  inner  (law  of  the  wall)  and  outer  (law  of  the  wake)  region 
defined  as  follows: 


Inner  Region: 


Re3/2  l 2  It  I 


(26) 


where 

*■  =  mixing  length  =  0.41  y  [1  -  exp  (~y+/A+)]  (27) 

A+  =  flat  plate  damping  length  constant  =  26.0 

+  V 

y  =  Reynolds  number  based  on  the  friction  velocity  -  -  (28) 

1/2 

u  =  friction  velocity  =  (x  /p)  (29) 

T  W 

Normally,  the  damping  length  constant  is  a  function  of  the  pressure 
gradient  at  the  wall,  (3p/3x)w.  Due  to  the  present  method  of  calculation, 
inclusion  of  this  term  presents  difficulties  that  are  too  time-consuming 
to  be  considered  for  this  report. 


Outer  Region: 
where 


Re  K  C  F„t„„ 
cp  KLEB 


(y)  f 


wake 


K 

C 

cp 

*wake 


Clauser  constant  -  0.0168 
1.60 


y  F  Re 
max  max 


-1/2 


0.25  y  u2  Re1/2/F 

max  DIF  max 


the  smaller 


F  and  u___  are  determined  from 

max  DIF 


(30) 
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F(y)  =  y  | t |  [1  -  exp  (-y+/A+)] 


(31) 


and 


UDIF  ~  °MAX  "  “MIN 


(32) 


FvTK-n(y)  =  Klebanoff  intermittancy  =  [1  +  5.5  (y  0.3/y  )6]-1  (33) 

MitD  IQ3X 


NUMERICAL  METHODS 

As  stated  earlier,  the  final  set  of  boundary-layer  equations  is  only 
weakly  elliptic.  The  normal  momentum  equation  can  then  be  uncoupled,  leaving 
the  conventional  boundary-layer  equations  which  are  solved  by  a  marching 
technique.  The  calculation  of  the  normal  pressure  gradient  becomes  an 
extra  step  in  the  strong-interaction  procedure,  lagging  one  cycle  in  the 
iteration.  The  other  major  modification  to  the  procedure  is  the  patching 
of  the  boundary-layer  and  potential-flow  solutions  at  the  boundary-layer 
edge  instead  of  matching  them  at  the  displacement  body  surface.  The 
numerical  and  computational  methods  required  to  institute  these  changes, 
including  smoothing  and  under-relaxation,  will  be  discussed  in  the  following 
sections. 


3.1  Boundary-Layer  and  Potential -Flow  Programs 

The  two  major  programs  employed  are  a  boundary-layer  program  and  an 
inviscid-flow  program.  The  boundary-layer  program  used  was  developed  by 
Hoffman  [ 8]  and  is  similar  in  method  to  that  presented  by  Cebeci  and 
Smith  [9].  It  is  a  general  finite-difference  solution  of  the  incompressible 
boundary-layer  equations  on  a  body  of  revolution  based  on  the  Keller  Box 
Method.  The  program  is  written  in  physical  variables,  necessitating  the 
input  of  a  starting  profile  as  a  singularity  is  encountered  at  the  body 
nose  in  this  formulation.  It  does,  however,  allow  for  calculations  to 
proceed  uninterrupted  into  the  wake,  a  requirement  in  this  elliptic 
problem  to  correctly  predict  the  flow  at  the  aft  end  of  the  body.  Programs 
using  the  Mangier  transformation  employ  an  integral  method  into  the  wake 
due  to  a  singularity  at  the  tail.  This  would  not  be  conducive  to  the 
present  investigation  because  discrete  values  of  the  flow  parameters  are 
required  throughout  the  boundary-layer  and  wake  flows  for  use  in  subsequent 
normal  pressure  calculations.  The  program  is  capable  of  treating  laminar 
or  turbulent  flow  up  to  a  separation  point  and  includes  transverse 
curvature  effects.  Major  modifications  to  the  program  included: 
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A.  Reformulation  of  the  boundary-layer  equations  in  terms  of  the  pressure 

and  addition  of  a  subroutine  that  calculates  the  pressure  variation 
across  the  boundary  layer  at  each  station.  The  pressure  gradient  is 
determined  by  applying  a  two-point-backwards  finite-difference  method 
to  the  normal  momentum  equation  and  using  values  of  u  and  V  saved  from 
the  previous  iteration. 

B.  Calculation  of  the  normal  component  of  the  velocity  at  each  station 
using  a  two-point-backwards,  finite-difference  method.  This  includes 

the  input  of  two  profiles  from  the  program  in  Mangler-Levy-Lees  variables. 

C.  Addition  of  a  subroutine  that  calculates  the  boundary-layer  thickness 
6  at  the  point  in  the  flow  where  u  *  0.995  ue.  Values  of  the  stream 
function  and  radius  are  then  determined  at  y  =  6.  These  values  define 
the  so-called  "displaced,  leaky  body,"  a  pseudo-body  consisting  of 

the  physical  body  with  the  boundary-layer  thickness  added  with  flow 
through  its  "walls"  from  boundary-layer  entrainment.  The  leaky  body  is 
used  as  input  to  the  potential  flow  program. 

D.  Addition  of  a  subroutine  that  changes  the  axial  step  size  from  a 
variable  to  a  constant  value  and  interpolates  for  ij^  and  at  the  new 
locations. 

The  inviscid-flow  program  was  also  developed  by  Hoffman  [5].  It 
involves  the  frozen  vorticity  equations  for  axisymmetric  flow  formulated 
using  a  spline-relaxation  technique.  The  Douglas-Neumann  program  for 
inviscid  flow  was  a  possible  option  and  in  fact  is  used  to  obtain  the 
initial  distribution.  Hoffman's  program  is  used  in  preference  for  later 
iterations  due  to  its  availability  and  the  ease  with  which  it  could  be 
modified  to  meet  the  needs  of  the  present  investigation.  Input  to  the 
program  includes  an  initial  profile  in  spline  variables  and  the  "leaky 
body"  information.  The  Cp  distribution  at  6  is  then  calculated  by  sweeping 
from  left  to  right  in  the  streamwise  direction  and  performing  line  relax¬ 
ation  at  each  station  until  a  converged  solution  is  obtained.  The  only 
modifications  to  the  program  involved  elimination  of  a  junction  region 
between  body  and  wake  calculations  that  was  made  irrelevant  by  the  cusped 
tail  on  the  F-57  body.  One  limitation  imposed  by  this  program  was  that 
the  input  data  must  be  in  constant  increments  in  the  axial  direction.  This 
was  accomplished  by  taking  the  variable  step  size  data,  fitting  it  with 
cubic  splines,  and  interpolating  for  new  values  with  a  constant  increment. 

Assorted  auxiliary  programs  are  used  between  runs  of  the  two  main 
programs  to  examine,  compare,  and  smooth  data.  Cubic  splines  are  employed 
in  accomplishing  the  latter  with  manual  control  over  the  choice  of  the 
spline  nodes. 
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3. 2  Solution  of  the  Normal  Momentum  Equation 

The  final  form  of  the  normal  momentum  equation  expressed  in  Eq.  (22) 

was 


1  3(v2) 

2  3y 


+  Re 


lE  = 

3y 


The  first  attempt  at  numerical  solution  of  this  equation  employed  the 
Keller  box  method,  the  basis  for  the  conventional  boundary  layer  portion 
of  the  program.  This  formulation  was  abandoned  without  a  test  run  when  an 
oscillation  appeared  in  the  calculation  of  v  which  employed  the  same 
method.  In  order  to  avoid  these  oscillations,  which  seem  to  be  inherent 
in  the  box  method,  the  derivative  of  the  normal  velocity  component  with 
respect  to  x  was  computed  using  a  second-order,  finite-difference  method 
with  a  variable  step  size.  The  general  form  of  this  equation  is 


i  +  1 


Ax 

f  1 

l-fi] 

— 

■sr +  2 

1  1 

K  ♦  i  - 'J 

Ax^^  +  Ax  ^ 


(34) 


where  f  is  the  variable  in  question,  i  -  1,  i,  and  i  +  1  are  three 
consecutive  stations,  and  Ax^  and  AX2  are  the  steps  in  x  between  the 
stations.  Centering  the  normal  momentum  equation  of  i  +  1  and  j  +  1/2, 
then  using  standard  central  differences  for  the  y-derivatives  gives  the 
following  result: 


(35) 


Equation  (35)  is  then  solved  downwards  from  the  boundary-layer  edge  employing 
a  given  edge  pressure  and  known  velocity  profiles.  The  equation's  final 
form,  with  the  i  +  1  subscript  understood,  is 
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Equation  (36)  is  allowed  to  lag  one  cycle  in  the  iterative  procedure. 

Values  of  u  and  v  are  saved  from  the  boundary-layer  run  and  are  used  with 
the  newly  calculated  edge  pressure  from  the  inviscid  program  to  calculate 
the  pressure  variation  through  the  boundary  layer  in  the  next  iteration. 

The  values  of  u  are  determined  as  one  of  the  variables  in  the  conventional 
boundary  layer  calculations.  The  method  for  calculating  v  will  be  considered 
next. 


3.3  Normal  Velocity  Calculations 

The  variation  of  the  normal  component  of  velocity  through  the  boundary 
layer  is  obtained  by  solving  for  v  from  the  continuity  equation.  An 
alternate  solution  involving  the  definition  of  the  Stokes  stream  function, 


v 


Ill 

r  3x 


(37) 


was  eventually  used  as  a  consistency  check,  but  the  continuity  formulation 
was  preferred  due  to  the  ease  with  which  it  could  be  extended  beyond  the 
boundary-layer  edge.  Beginning  with  the  continuity  equation,  we  have 


~  (ru)  +  ■—  (rv)  =  0  .  (38) 

3x  3y 

As  stated  earlier,  the  use  of  the  Keller  box  method  in  solving  this  equation 
resulted  in  oscillations  in  the  profiles  from  station  to  station.  The 
problem  was  resolved,  as  in  the  pressure  calculation,  by  using  backward 
differences  for  the  derivatives  in  the  x-direction  and  then  centering 
along  a  vertical  element  of  the  boundary-layer  grid.  The  resulting  equation, 
solved  for  (rv)j+1  at  station  i  +  1,  is 


(rv) 


j+1 


3__ 

3x 


(rU)J+i 


(39) 


Equation  (39)  is  solved  from  the  body  surface,  where  v  =  0,  upwards  to 
the  boundary-layer  edge.  J 


3.4  Determination  of  the  Boundary-Layer  Edge  and  Associated  Values 

The  boundary-layer  and  inviscid-f low  solutions  are  patched  together 
at  the  boundary-layer  edge  in  the  present  calculation.  Therefore  the 
location  of  the  boundary-layer  edge  and  the  values  of  several  parameters 
at  the  edge  must  be  determined.  The  values  of  the  radius  and  Stokes 
stream  function  at  6  are  required  as  input  to  the  potential-flow  program, 
v  is  also  needed  to  calculate  ue  in  the  next  iteration  of  the  boundary- 
layer  program.  It  is  desirable  to  have  these  distributions  as  smooth  as 
possible  when  first  calculated  to  reduce  the  amount  of  manual  smoothing 
required.  Several  methods  of  determining  the  boundary  layer  were  explored 
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and  abandoned  due  to  unsatisfactory  results.  Each  will  be  briefly 
described  with  the  results,  which  led  to  the  best  choice,  shown  in  terms 
of  and  r^  in  Figure  3. 

The  initial  procedure  was  simply  to  allow  the  parameters  to  be  defined 
at  the  last  point  in  the  boundary-layer  calculation  at  each  station.  Due 
to  the  nature  of  these  calculations,  specifically  the  choppy  growth  pattern 
of  the  boundary-layer  "edge,"  the  results  are  very  ragged  and  therefore 
highly  unsatisfactory. 

The  next  method  involved  evaluation  of  6  at  the  y  value  where 
u  -  0.995  u  .  The  interpolation  was  accomplished  with  cubic  splines. 

This  approximate  boundary-layer  edge  was  then  used  to  determine  r^, 
and  vg  by  a  similar  interpolation.  This  method  resulted  in  a  relatively 
smooth  distribution  except  where  oscillations  in  u  occurred.  Note  that 
these  oscillations  occur  in  the  tail  region  of  the  body.  The  fluctuations 
again  make  this  method  undesirable. 

The  final  method  was  to  input  a  fixed,  smooth  curve  as  the  boundary- 
layer  edge.  The  curve  was  initially  determined  by  taking  the  output  from 
the  second  method  and  manually  smoothing  it  with  a  cubic-spline  routine. 
Figure  4  shows  the  boundary-layer  edge  determined  in  this  way.  It  is 
input  into  every  iteration  of  the  boundary-layer  program  and  the  values 
of  r^,  ipfl,  and  ve  are  calculated  at  the  corresponding  point  in  the  flow. 

The  results  were  generally  very  good.  Examining  Figure  5,  it  can  be 
seen  that  each  distribution  is  at  least  as  smooth  as  the  curve  used  for 
interpolation.  This  method  allows  for  ease  in  patching  the  potential 
and  boundary-layer  solutions  at  the  boundary-layer  edge.  It  has  been 
suggested  [10] , however ,  that  the  approximated  boundary-layer  edge  may 
have  to  be  revised  in  later  iterations. 


3.5  Calculations  Beyond  the  Boundary-Layer  Edge 

The  box  method  employed  in  the  boundary-layer  program  requires  the 
extension  of  the  boundary-layer  edge  when  the  value  of  the  shear  stress 
increases  over  some  cutoff  value.  The  nature  of  the  normal  velocity  and 
pressure  calculations  requires  that  there  be  grid  points  and  associated 
values  for  two  stations  immediately  behind  the  station  in  question.  Thus, 
when  the  boundary-layer  edge  is  extended  at  station  i  +  1,  v  and  p  must 
be  obtained  beyond  the  boundary  layer  at  stations  i  and  i  -  1. 

First  of  all,  the  simplifying  assumption  that 


u  =  u  for  y  >  6 
e  — 


(40) 


is  made.  Then  v  can  be  determined  from  the  continuity  equation  where 
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W  (rv)  -  *  h  (->  -  -  h  (rue)  • 

Integrating  Equation  (41)  from  y  to  y,+1  and  using  the  definition  of 
r(x),  we  obtain  the  final  form  or  the  equation: 


(rv)i,j+l  »  (rv)i 


6(xi)  fc+1  +  y1 

.j  -  a(xt)  +  ~ —  j ^  Ayi 

/Re  3 


where 


a(xi}  *  ¥x  (ue  ro(x))  ’ 


e(xi>  ■  (Ue  COS*(x)) 


ayJ  ■  yJ+l  •  yJ  • 

v  can  be  determined  directly  from  Equation  (42)  as  the  value  of  r  is 

known. 

Since  we  are  above  the  boundary-layer  edge  and  in  potential  flow, 
Bernoulli's  equation  holds.  Written  dimensionally,  it  is 

*  1  *  *2 

p  +  -^p  V  ■  constant  ■  c^  (43) 


If  we  define  dimensionless  variables  as  follows: 


*  *  *2 
P  ■  P  U  p  , 


*  * 

V  -  U  V  , 


then  Equation  (43)  in  dimensionless  form  is 


.  1  ..2  C1 

P  +  2  V  =  -T-^2  m  C2  * 

P 


constant  c^  can  be  evaluated  at  the  boundary  layer  edge  as  follows: 

c2  "  pe  +  1  (ue2  +  ve2)  •  <4' 
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« 
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Writing  Equation  (45)  for  a  point  above  the  boundary-layer  edge  and 
substituting  Equation  (46)  for  c ^  we  have 

12  2  12  2 
p+-r(u  +v)=p  +-r(u  +v  )  .  (47) 

r  2  e  re  2  e  e 

Solving  for  p  and  rewriting  in  terms  of  the  scaled  velocity,  v,  we  have 
finally 


Pi,j+1 


Vi,j+1  > 


(48) 


The  value  of  the  pressure  above  the  boundary  layer  can  now  be  calculated 
employing  the  value  of  V.  .  -  determined  from  Equation  (42) . 

*  J  j  ' 


3.6  Calculation  Loops 

With  the  equations  and  boundary  conditions  as  specified  in  section  2, 
the  boundary  layer  calculations  proceed  as  follows: 

1.  Given  Cp5  at  station  x^,  the  normal  momentum  equation  is  solved  down¬ 
wards  to  the  body  using  tt e  values  of  u(x,y)  and  v(x,y)  from  the 
previous  iteration.  The  first  iteration  is  obtained  by  assuming  that 
3p/3y  =  0. 

2.  Holding  3p/3y  fixed,  the  conventional  boundary  layer  equations  are 
solved  from  which  values  of  and  r^  are  determined. 

3.  The  normal  velocity  profile,  v(x,y)  is  calculated  and  stored  along 
with  the  streamwise  velocity  profile  for  use  in  the  next  iteration. 

4.  Steps  (1  -  3)  are  repeated,  marching  downstream  until  the  desired 
flow  field  is  obtained. 

Both  the  normal  and  streamwise  components  of  velocity  must  be  stored  in 
two-dimensional  arrays.  The  necessary  storage  can  be  minimized  by  including 
the  pressure  calculation  only  where  the  normal  pressure  gradient  is 
significant.  This  portion  of  the  flow  field  is  from  about  70%  of  the 
body  length  onwards,  as  can  be  seen  in  Figure  8.  Present  calculations 
were  carried  out  beginning  at  about  26%  of  the  body  length  as  a  test  of 
the  accuracy  of  the  method.  There  should  be  very  little  variation  between 
the  results  using  the  present  method  and  those  using  the  conventional 
method  (which  neglects  normal  pressure  gradient  variations)  for  this  early 
region. 

The  calculation  loop  for  the  entire  modified  strong-interaction  pro¬ 
cedure  is: 

a.  Determine  an  initial  guess  for  C^(x),  including  into  the  wake  region. 
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b.  Compute  the  boundary-layer  flow  with  3p/3y  *  0,  giving  6(x),  r^x), 
and  ip  p (x) •  This  effectively  defines  the  "leaky  body."  Values  of 
v(x,y)  and  u(x,y)  are  determined  and  stored  for  the  next  iteration. 

c.  Calculate  Cpg(x)  from  the  "leaky  body"  inviscid  solution. 

d.  Compute  the  boundary  layer  flow,  calculating  the  normal  pressure 
variation  at  each  station  using  .CP6  (x)  determined  in  step  (c)  and 
the  velocity  profiles  from  previous  iterations.  With  the  values  of 
6(x),  determined  in  step  (b),  calculate  new  values  of  t|».(x)  and 

*«<*). 

e.  Repeat  steps  (c)  and  (d)  until  convergence  is  achieved. 

A  flow  chart  of  the  procedure  is  presented  in  Figure  6. 


3.7  Smoothing  and  Under-Relaxation 

The  numerical  representation  of  distinct  values  of  C  and  ^  requires 
some  care  as  the  boundary-layer  and  invisc id-flow  solutions  depend  on  the 
streamwise  derivatives  of  these  quantities.  Oscillations  in  either  would 
be  amplified  and  propagated  through  succeeding  calculations.  It  should  be 
noted,  however,  that  the  inviscid-flow  formulation  being  used  in  this 
Investigation  is  less  sensitive  to  variations  in  ip  than  conventional  pro¬ 
grams.  Some  smoothing  of  these  and  other  calculated  values  seems  to  be  in 
order. 

Significant  effort  was  expended  to  obtain  relatively  smooth  primary 
output  as  described  in  Section  3.4  for  the  boundary-layer  edge  calculation. 
The  amount  of  extra  smoothing  between  programs  has  so  far  been  restricted 
to  the  values  of  Cp^  from  the  potential-flow  program  because  its  x-wise 
stations  do  not  coincide  with  those  of  the  boundary-layer  program. 

A  sweep-to-sweep  under-relaxation  factor  of  0.3  applied  to  v  and 
C  ,  has  been  suggested  by  Mahgoub  and  Bradshaw  [10J.  As  will  be  shown  in 
tne  following  sections,  a  value  of  0.1  seems  to  be  more  appropriate  for 
the  first  few  iterations.  Convergence  should  be  at  least  monotonic  and 
could  possibly  be  improved  by  increasing  the  under-relaxation  factor  for 
later  iterations  [10]. 


RESULTS  AND  DISCUSSION 


4.1  Verification  of  Computational  Methods 

Validity  checks  were  made  on  the  newly  developed  components  of  the 
strong-interaction  procedure.  Both  main  programs  delivered  accurate 
results  when  run  as  conventional  programs,  that  is,  with  all  modifications 
turned  off.  Then  the  two  major  additions,  the  calculation  of  v  and  p 
through  the  boundary  layer,  were  checked  for  consistency. 


i 
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Calculations  of  v  using  the  continuity  equation  were  checked  against 
calculations  of  v  using  the  definition  of  the  stream  function, 


v 


1 

r  3x 


The  difference  between  values  obtained  from  the  two  methods  was  less  than 
0.2%.  Further,  upon  comparing  results  of  the  present  calculation  with  the 
theoretical  values  obtained  by  Patel  and  Lee  [3],  it  was  found  that  both 
methods  provide  similar  profiles,  but  they  are  much  too  full.  The  ratio 
of  v@  to  in  both  cases  is  approximately 

v 

theor 

as  can  be  seen  in  Figure  7.  Each  calculation  was  made  with  3p/3y  =  0. 
Better  normal  velocity  profiles  should  be  obtained  in  later  iterations 
of  the  present  procedure  when  the  pressure  is  varied  through  the  boundary 
layer . 


The  method  of  determining  the  normal  pressure  gradient  was  also 
checked.  The  profiles  obtained  corresponded  to  experimental  values  up 
to  about  70%  of  the  body  length.  At  this  point,  the  values  of  v-  and 
3v/3x  employed  in  the  normal  momentum  equation  both  become  much  larger 
than  experimental  measurements  indicate.  Consequently,  the  pressure 
calculations  from  this  point  on  are  inaccurate  and  resulting  profiles 
do  not  agree  with  experimental  values.  Resolution  of  this  difficulty 
is  the  major  stumbling  block  left  in  the  investigation  and  is  the  target 
of  present  attempts  at  under-relaxation  and  smoothing. 

The  ability  of  the  patched  viscous-inviscid  solution  to  correctly 
predict  the  edge  pressure  was  examined  next.  The  experimental  values  of 
Cp  at  the  wall  were  used  as  input  data  to  the  boundary-layer  program. 
Values  of  Cp.  corresponding  to  experimental  measurements  were  the  expected 
output  from  the  potential-flow  program.  Results  are  shown  in  Figure  8. 

The  pressure  distribution  obtained  follows  the  experimental  data  very 
closely,  even  modeling  the  dip  in  the  tail  region.  The  calculated  values 
are  slightly  lower  than  experiment  in  this  region,  probably  due  to  the 
boundary-layer  edge  approximation  employed. 

One  result  that  appeared  to  be  incorrect  was  the  dip  that  occurs  at 
about  80-95%  of  the  body  length  in  the  curve  of  tjjg  vs.  xg  as  shown  in 
Figure  5.  The  curve  follows  a  smooth,  increasing  trend  everywhere  except 
in  this  area  and  the  fluctuation  is  not  reflected  in  any  of  the  other 
boundary-layer  parameters.  The  validity  of  this  dip  was  checked  in  the 
following  manner. 

The  rate  of  change  of  the  stream  function  with  respect  to  x  at  the 
boundary  layer  edge  can  be  written  as 
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'd£ 

dx 


6 


From  the  definition  of  the  stream  function  we  have 


v 

e 


(49) 


(50) 


Substituting  Equation  (50)  into  Equation  (49)  we  obtain 


dx 


i*l] 

dxl 

'  t 


(51) 


Equation  (51)  allows  a  direct  examination  of  the  slope  of  the  curve  of 
and  insight  into  the  validity  of  the  dip.  All  values  on  the  right- 
hand-side  of  the  equation  are  output  from  the  boundary-layer  program. 
Magnitudes  of  the  terms  can  be  obtained  from  Figures  4,  5,  and  9.  It  can 
be  seen  that  the  normal  component  of  the  velocity  v  increases  rapidly  to 
a  peak  at  approximately  90%  of  the  body  length,  whife  at  the  same  point 
r.  decreases  to  its  lowest  value  for  the  rear  portion  of  the  body.  From 
tne  plot  of  6,  tne  boundary-layer  thickness,  we  can  see  that  (dy/dx)^  is 
constant  in  the  vicinity  of  the  tail  of  the  body.  The  value  of  u  decreases 
as  the  flow  approaches  the  aft  end  of  the  body.  The  combined  effict  points 
to  a  decreasing  d^/dx  over  the  range  from  80%  to  90%  and  gives  a  plausible 

explanation  to  the  dip  in  ip  _  that  occurs  in  this  area. 

0 

Overall,  the  findings  of  these  tests  have  strengthened  the  basis  for 
this  investigation.  They  have  proven  that  the  pressure  distribution  can 
be  calculated  accurately  at  the  boundary-layer  edge,  allowing  for  the 
patching  of  the  viscous-inviscid  flow  solutions  there.  They  have  shown 
that  the  normal  velocity  and  pressure  variation  through  the  boundary 
layer  can  be  calculated.  Most  importantly,  they  point  the  way  for  the 
additional  numerical  experimentation  required  to  overcome  current  dif¬ 
ficulties. 


4.2  Convergence 

A  converged  pressure  distribution  is  the  ultimate  goal  of  the  investi¬ 
gation.  To  achieve  this,  several  problems  must  be  resolved.  By  examining 
Figure  10,  it  is  seen  that  the  C  distribution  used  as  a  "first  guess" 

wall 

for  the  boundary-layer  program  differs  greatly  from  the  experimental 
values.  This  results  in  the  profiles  of  v  being  too  full  at  the  tail,  as 
shown  in  Figure  7,  and  the  large  over-correction  in  C  obtained  from  the 
inviscid  flow  program  as  seen  in  Figure  10.  ^5 
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Two  basic  difficulties  in  the  problem  are  the  extreme  sensitivity  of 
the  calculations  at  the  tail  and  the  oscillations  in  u  in  this  region. 

The  flow  field  changes  very  rapidly  due  to  the  large,  streamline  curvature 
at  the  tail.  Small  errors  in  the  calculations  tend  to  be  magnified  thus 
producing  poor  convergence  characteristics.  To  minimize  the  change, 
under-relaxation  factors  on  the  order  of  0.1  seem  to  be  called  for  in 
the  calculation  of  the  normal-velocity  and  normal-pressure  distributions 
in  the  first  few  iterations.  The  oscillations  in  u  at  the  tail  of  the 
body  are  inherent  in  the  present  physical  variable  finite  difference 
formulation  of  the  boundary-layer  program.  Hopefully,  they  will  be 
smoothed  out  as  the  iteration  cycle  proceeds  and  the  flow  field  model 
becomes  more  precise. 

Stated  simply,  the  ultimate  convergence  of  the  present  investigation 
depends  heavily  on  the  successful  application  of  under-relaxation  techniques. 


4.3  Conclusions 


A  method  for  including  the  normal  pressure  gradient  in  the  axisymmetric , 
strong-interaction  calculations  has  been  presented.  The  determination  of 
the  pressure  variation  through  the  boundary  layer  via  solution  of  the 
normal  momentum  equation  and  the  patching  of  the  boundary-layer  and 
potential  flow  solutions  at  the  boundary-layer  edge  are  the  major  innovations. 
The  procedure  was  applied  to  the  F-57  body,  a  low-drag  body  of  revolution 
with  a  cusped  tail.  A  converged  solution  has  not  been  achieved  as  of  this 
writing.  Sensitivity  of  the  calculations  along  with  the  large  error  in 
the  initial  pressure  profile  have  limited  success  in  this  direction.  Current 
research  into  the  problem  is  centered  around  under-relaxation  techniques. 
Proper  application  of  these  techniques  should  result  in  a  converged  solution. 

The  viability  of  the  method  has  been  at  least  partially  established 
by  several  of  the  tests  conducted,  thus  creating  a  basis  for  further 
investigation.  Continued  research  should  produce  a  more  accurate  means 
of  modeling  the  flow  field  at  the  aft  end  of  a  body  of  revolution. 
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Figure  4.  Boundary-Layer  Thickness,  F-57  Body. 
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Figure  6.  Interaction  Scheme  for  the  Modified  Strong-Interaction 
Method . 
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Figure  10,  Comparison  of  C„  on  the  Body  and  at  the  Boundary— Laver 
Edge,  F-57  Body, 
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